# TODO: Add comment
# 
# Author: broken
###############################################################################
#import library
library(base)
library(raster)
library(gstat)


#path

path= "C:/Users/forever/Desktop/MYD01/"

year_array=c(2010,2011,2012,2013,2014)
for (j in 1:5) {
	inputPath = paste(path,year_array[j],sep="")
	#pm file
	fileArray = list.files(path = inputPath,pattern="10km _og.tif$",full.names = FALSE,recursive=TRUE)
	totalFile=length(fileArray)
	
	## #co file
	## demFile = "E:/FIMO/DATA/PM/DEM/dem.tif"
	## nlFile = paste("E:/FIMO/DATA/PM/NL",year_array[j],"nl_10km.tif",sep="/")
	## lcFile = paste("E:/FIMO/DATA/PM/lC",year_array[j],"B1_10km.tif",sep="/")

	## #array
	## demCorArray= nlCorArray= lcCorArray= c(1:totalFile)
	## 
	## sCor=sRMSE=sRE=c(1:totalFile)
	## 
	## oCor=oRMSE=oRE=c(1:totalFile)
	
	uCor=uRMSE=uRE=numberPixel=c(1:totalFile)
	
	for(i in 1:totalFile){
		
		fileName=fileArray[i]
		
		#PM values
		pmFile=paste(inputPath,fileName,sep="/")
		pmRaster=raster(pmFile)
		pm=values(pmRaster)
		corxy=coordinates(pmRaster)
		x=corxy[,'x']
		y=corxy[,'y']
		
		
		## #dem values
		## demRaster=raster(demFile)
		## dem=values(demRaster)
		## 
		## #night light values
		## nlRaster=raster(nlFile)
		## nl=values(nlRaster)
		## 
		## #landcover values	
		## lcRaster=raster(lcFile)
		## lc=values(lcRaster)
		## 
		#cell id
		totalCell=length(pmRaster)
		cell = c(1:totalCell)
		
		#data frame
		#table=data.frame(cell,x,y,pm,dem,nl,lc)
		table=data.frame(cell,x,y,pm)
		newTable=table
		testTable=subset(table,pm<0)
		trainTable=subset(table,pm>=0)

		## #caculate corelation between pm,dem,nl,lc
		## corelation=cor(data.frame(trainTable$pm,trainTable$dem,trainTable$nl,trainTable$lc))
		## demCorArray[i]=corelation[1,2]
		## nlCorArray[i]=corelation[1,3]
		## lcCorArray[i]=corelation[1,4]
		
		#caculate variogram
		empiVario=variogram(pm~1,locations=~x+y,data=trainTable)
		
		#sph fit
		
		sill=min(empiVario$gamma)
		#sphModel=vgm(psill=sill,model="Sph",nugget=0,range=min(empiVario$dist))
		sphModel=vgm(model="Sph",nugget=0,range=1)		
		sphFit=fit.variogram(empiVario,sphModel)
		
		#plot sph variogram to jpeg
		## outFile=paste(inputPath,"/",substr(fileName,1,40),"_vario.jpg",sep="")
		## jpeg(file=outFile)
		## plot(empiVario,model=sphFit)
		## dev.off()
		## 
		## #-----------------------SIMPLE KRIGING-----------------------------------
		## 
		## beta_mean=mean(trainTable$pm)
		## simple_result=krige(id="pm",formula=pm~1,data=trainTable,newdata=newTable,model=sphFit,locations=~x+y,beta=beta_mean)
		## 
		## #edit tiff
		## simplePMRaster=pmRaster
		## simplePMValue=simple_result[,3]
		## simplePMRaster[1:totalCell]=simplePMValue
		## 
		## #save sk result to tiff
		## outFile=paste(inputPath,"/",substr(fileName,1,40),"_sk.tif",sep="")
		## writeRaster(simplePMRaster,filename=outFile,format="GTiff")
		## 
		## #cross-validation
		## simple_result_3=krige.cv(pm~1,trainTable,sphFit,locations=~x+y,nfold=3,beta=beta_mean)
		## 
		## 
		## #simple statis
		## simple_cor=cor(simple_result_3$var1.pred,simple_result_3$observed)
		## simple_rmse=sqrt(sum((simple_result_3$residual)^2)/nrow(simple_result_3))
		## simple_re=sum(abs(simple_result_3$residual)/simple_result_3$observed)/nrow(simple_result_3)
		## 
		## sCor[i]=simple_cor*simple_cor
		## sRMSE[i]=simple_rmse
		## sRE[i]=simple_re*100
		## 
		## print("simple completed")
		## 
		## 
		## #-----------------------ODINARY KRIGING-----------------------------------
		## 
		## odinary_result=krige(id="pm",formula=pm~1,data=trainTable,newdata=newTable,model=sphFit,locations=~x+y)
		## 
		## #edit tiff
		## odinaryPMRaster=pmRaster
		## odinaryPMValue=odinary_result[,3]
		## odinaryPMRaster[1:totalCell]=odinaryPMValue
		## 
		## #save ok result to tiff
		## outFile=paste(inputPath,"/",substr(fileName,1,40),"_ok.tif",sep="")
		## writeRaster(odinaryPMRaster,filename=outFile,format="GTiff")
		## 
		## #cross-validation
		## odinary_result_3=krige.cv(pm~1,trainTable,sphFit,locations=~x+y,nfold=3)
		## 
		## 
		## #Odinary statis
		## odinary_cor=cor(odinary_result_3$var1.pred,odinary_result_3$observed)
		## odinary_rmse=sqrt(sum((odinary_result_3$residual)^2)/nrow(odinary_result_3))
		## odinary_re=sum(abs(odinary_result_3$residual)/odinary_result_3$observed)/nrow(odinary_result_3)
		## 
		## oCor[i]=odinary_cor*odinary_cor
		## oRMSE[i]=odinary_rmse
		## oRE[i]=odinary_re*100
		## 
		## print("Odinary completed")
		
		
		#-----------------------UNIVERSAL KRIGING-----------------------------------
		
		universal_result=krige(id="pm",formula=pm~x+y,data=trainTable,newdata=newTable,model=sphFit,locations=~x+y)
		
			
		#edit tiff
		universalPMRaster=pmRaster
		universalPMValue=universal_result[,3]
		universalPMRaster[1:totalCell]=universalPMValue
		
		#save uk result to tiff
		outFile=paste(inputPath,"/",substr(fileName,1,40),"_uk.tif",sep="")
		writeRaster(universalPMRaster,filename=outFile,format="GTiff")
		
		#edit error tiff
		errorPMRaster=pmRaster
		errorPMValue=universal_result[,4]
		errorPMRaster[1:totalCell]=errorPMValue
		
		#save uk result to tiff
		outFile=paste(inputPath,"/",substr(fileName,1,40),"_error.tif",sep="")
		writeRaster(errorPMRaster,filename=outFile,format="GTiff")
		
		#cross-validation
		universal_result_3=krige.cv(pm~x+y,trainTable,sphFit,locations=~x+y,nfold=3)
		
		
		#Universal statis
		universal_cor=cor(universal_result_3$var1.pred,universal_result_3$observed)
		universal_rmse=sqrt(sum((universal_result_3$residual)^2)/nrow(universal_result_3))
		universal_re=sum(abs(universal_result_3$residual)/universal_result_3$observed)/nrow(universal_result_3)
		numberPixel[i]=nrow(trainTable)
		uCor[i]=universal_cor*universal_cor
		uRMSE[i]=universal_rmse
		uRE[i]=universal_re*100
		
		print("Universal completed")
		
		
		## #plot all to jpeg
		## pmOrigi=pmRaster
		## pmOrigi[pmOrigi<0]=NA
		## outFile=paste(inputPath,"/",substr(fileName,1,40),"_all_compare.jpg",sep="")
		## jpeg(file=outFile)
		## par(mfrow=c(2,2))
		## plot(pmOrigi,main="Original")
		## plot(simplePMRaster,main="SK")
		## plot(odinaryPMRaster,main="OK")
		## plot(universalPMRaster,main="UK")
		## dev.off()
		## 
		
		str=paste("Complete",fileArray[i],sep=" ")
		print(str)
	}
	

#export statis result
	result_dataframe=data.frame(fileArray,numberPixel,uCor,uRMSE,uRE)
	outFile=paste(inputPath,"/","cross_result.csv",sep="")
	write.csv(result_dataframe,file=outFile)
	print(paste("Finish",year_array[j]))
	
	

}






